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Preface 


About 6 years ago, one of the authors (PFB) was faced with 
the requirement of performing spectral calculations for ill-behaved 
waveforms that cannot be easily treated using the Fast Fourier 
transform. Examples of this type of waveform include a unit step 
function and a semi-infinite sine wave. 


He had heard anecdotal stories of some researchers modifying 
this type of waveform by multiplying it by an attenuating 
exponential function exp(—o,f) to artificially dampen the 


waveform so that a numerical Fourier transform could be taken. 
Then, after manipulating the computed spectral response (say, by 
applying a known transfer function, or by time shifting the 
spectrum), an inverse Fourier transform is performed. The resulting 
transient waveform is then multiplied by exp(+o,t) to remove the 


effects of the artificially induced damping function. 


In thinking about this procedure, he realized that this is 
nothing more than a description of the Lap/ace transform of the 
function f(/), which yields the spectrum not at a frequency f= @/27, 
but at a complex frequency s = o, + j@. Bertholet went on to 


develop a simple direct and inverse Laplace transform routine that 
uses the FFT algorithm to perform the necessary calculations. This 
concept was documented internally for his laboratory in French, but 
it was never published in the open literature. 


Two papers in the electrical power community in 2004 and 
2007 have described this computational procedure and have called 
it the numerical Laplace transform (NLT). These references have 
provided a historical review of the past work leading to the 
realization of this approach, and have given several illustrations of 
the use of this method for power system analysis. This method is 
actually based on an earlier paper in 1968, which refers to the 


method as the modified Fourier transform (MFT). 


These earlier references provide a review of the theoretical 
basis for the method of computing the Laplace spectrum of f(/) and 
in computing its inverse. However, the details of how the Laplace 
spectrum is to be used for practical problems are not discussed 
explicitly in this earlier work. This monograph serves to fill this 
void. 


This monograph reviews the use of the Laplace transform as 
implemented using the fast Fourier transform. While this method is 
used in the electrical power community, it is not widely used in the 
electromagnetic compatibility (EMC) area. The goal in developing 
this monograph, therefore, is to bring this computational method to 
the attention of the workers in the EMC community by providing 
several examples and comments on its use for practical problems. 


We ate indebted to the Swiss Federal Office for Civil 
Protection of the Federal Department of Defense, Civil Protection 
and Sports for their support of the investigations leading to this 
document. In particular, we would like to acknowledge discussions 
and comments from Dr. Alain Jaquier and Mr. Markus Nyffeler, 
both from the Swiss armasuisse organization. 
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Chapter 1. Introduction 


The numerical Fourier transform and its inverse transform are 
powerful analysis tools often used in the electromagnetics (EM) 
area. The development of the fast Fourier transform (FFT) by 
Gauss citca 1805 [1] and the subsequent machine implementation 
by Cooley and Tukey some 165 years later [2], have contributed 
greatly to our ability to pass from the time domain to the frequency 
domain, and back, with very large data records. 


There are limitations, however, to Fourier analysis methods, 
which ate usually encountered by anyone who tries to use them. A 
sufficient condition for a Fourier transform of a time-domain 
function f(/) to exist, is that 4 must be absolutely integrable [3]: 


ioe) 


| |f@|dt<@. (1) 


—oO 


This condition is satisfied by all practical transient signals that have 
a finite duration’. However, some idealized waveforms, like the step 
function and the semi-infinite sinusoidal waveform, pose difficulties 
for Fourier transformation. 


In engineering applications, the Fourier transform of a 
waveform can be computed either by a direct evaluation of the 
Fourier integral, or by the use of the discrete Fourier transform 
(DFT). In this latter approach, the transient waveform and the 
resulting spectral response are sampled at discrete time and 
frequency points. The FFT is a special case of the DFT, in which 
the number of points in the transient waveform, N, is constrained 


to be equal to N = 2”, where m is an integer. 


1 As noted in [3], (1) is a sufficient, but not necessary condition for the existence of the 
Fourier transform of f(/. Papoulis give as an example the function f/ = sin(@A)/4, which 
does not satisfy (1), but which does have a Fourier transform. 


In practical cases where numerical evaluations are used, the 
time window of the transient signal must be truncated. In addition, 
the discrete sampling of f# implies that the Fourier spectrum is 
band limited. Thus, the numerical transform is only an 
approximation to the actual Fourier transform of a waveform. If 
the sampling of f(t) and the time window in which the waveform is 
defined are not properly chosen, the computed spectrum can be 
incorrect. 


These observations have led to two practical guidelines for 
performing numerical Fourier transforms, which are 


e Choose the time window in which the waveform f(/) 
is defined so that the waveform is zero at each end of 
the window, and 


e Choose the number of points in the discrete sampling 
of f(4) so that the rate of change of the fastest part of 
the waveform is adequately represented. 


While these are seemingly rather straightforward requirements, 
they occasionally pose problems for computing the transform of 
certain types of functions. Several of these functions are shown in 
Figure 1. This figure shows a commonly used model for the return 
stroke current in a lightning channel and the resulting E-field [4], a 
step function and a semi-infinite sinusoidal signal. Each waveform 
is started at a different time for clarity, and the amplitudes are in 
arbitrary units. 


In addition to these signals, waveforms that are distributions, 


such as the Dirac delta function 8(/ or f4 = 1/ vi , pose problems 
with the FFT because it is impossible to represent them numerically 
as discretely sampled functions’. 


2 Of course, these functions also can be handled analytically. 


Lightning channel Lightning-produced 
current 


f(t) 


Sinusoidal signal 


0 20 40 60 80 100 
Time (us) 


Figure 1. Illustration of several simple waveforms that pose 
difficulties in computing the Fourier transform using the FFT. 


To perform spectral calculations of ill-behaved waveforms like 
those in Figure 1, it is possible to use the same method that is 
employed with analytically passing from the Fourier transform to 
the Laplace transform. That is, we can modify the initial waveform 


by an attenuating exponential function exp(—o,t) to artificially 


dampen the waveform in the time window so that a Fourier 
transform can be taken. Then, upon the manipulation of the 
computed spectral response, an inverse Fourier transform can be 
performed. The resulting transient waveform is then multiplied by 


exp(+o,, t) to remove the effects of the initial artificial damping. 


This procedure is essentially a description of the Laplace 
transform of the function f(4, which yields the spectrum, not at a 
real frequency f= @/27, but at a complex frequency s = o, + ja. 
With this method, it is possible to develop a simple direct and 
inverse Laplace transform routine that uses the FFT algorithm to 
perform the necessary calculations. 


Recently, two papers in the electric power community have 
described this computational procedure, and have referred to it as 
the numerical Laplace transform (NLT) [5, 6]. These references have 
provided a historical review of the past work leading to the 
realization of this approach and have given several illustrations of 


the use of this method for power system analysis. This method has 
its roots in an earlier paper [7], which calls this method the modified 
Fourier transform (MFT). 


These earlier references provide a good review of the 
theoretical basis for the method of computing the Laplace spectrum 
of {4 and in computing its inverse. However, the details of how the 
computed spectrum should be used for practical problems are not 
discussed explicitly. For example, if a signal f(/) is passed through a 
linear system (say a filter) its Laplace spectrum must be multiplied 
by the Laplace spectrum of the system’s transfer function evaluated 
at the complex frequency (s) rather than at the real frequency (/) in 
the Fourier domain. While this does not pose a problem for transfer 
functions that are known analytically, care must be used if the 
transfer function is measured in the Fourier domain. Such transfer 
functions in the Fourier domain must be analytically continued into 
the complex frequency domain to be used with the Laplace 
transform. 


While the Laplace transform method of [5, 6] appears to be 
well-known to the power community, a poll of investigators in the 
EMC area suggests that this method is not commonly used by 
them, and that many people are unaware of its benefits. 


In this monograph, we will refer to this numerical procedure as 
the fast Laplace transform (FLT) to be consistent with the notion of 
the Laplace transform as commonly used in the electrical 
engineering literature. We describe this method as “fast” because 
the FFT algorithm is used in the process. 


Chapter 2. Theoretical 
Development 


2.1. The Fourier and Laplace Transforms 


For a transient waveform /(/ that is zero for ¢ < 0, its Fourier 
transform F(@) is given by [3] as 


F(@)= | f(Derdt . (2) 


As noted earlier, a sufficient (but not necessary) requirement for the 
existence of F(@) is given by (1). The inverse Fourier transform, 
which provides the transient waveform from the spectrum, is given 


by 


ro) == J Floyeao @) 


For waveforms that do not meet the requirement of (1), yet 
which are exponentially bounded, a spectrum can be computed by 
multiplying f/) by an exponential damping function to create a 
modified waveform f,,(t;O,) as 


fiGo,)=fOer* . (4) 


In this manner, the spectrum of this modified waveform becomes 


F.(@;0,) = i fe )e dt 
: (5) 
=| f.Go,)e?*at , 


where f,, is a function of both time and the damping constant Oo. In 
this expression, we note that the modified waveform /,, and the 


spectrum F,,, are related by the Fourier transform of (2). 


The inverse Fourier transform of the spectrum F,, provides 


the modified waveform f,, as 


1 t jot 
f,(to,) = 55 | Fal@so,)e do (6) 


and the transient function f(/) may be reconstructed by multiplying 
Inby 


f(=e"" us | F (a;0,)e"do . (7) 
Zit 


With s=o +jo@ and F(@;0,)=F(s), (5) and (7) form the 


well known Laplace transform pair” 


F(s)= { fe“dt , and 6) 
0 aaj | Povea 0) 


Note that (8) is the so-called one-sided, or unilateral, Laplace 
transform, and arises due to the fact that (4 for ¢< 0 is assumed to 
be zero. A more general expression for this transform exists in the 
form of the two-sided, or bilateral Laplace transform, as discussed 
in [8]. In this monograph, however, we will use only the one-sided 


3 The reader is cautioned that there is a possibility of confusion here with regard to the 
notation for the Fourier and Laplace transforms of the function f(/). In this monograph, 
and in the engineering literature, the function F(@) or F() (where @ = 27/) is commonly 
used to represent the Fourier transform of (2), and F(s) (where s = o,+j@) denotes the 


Laplace transform of Eq(8). 


transform, as all practical signals of interest are zero prior to some 
“turn-on” time, which we take to be at ¢= 0. 


The inverse Fourier transform of f# in (3) and the inverse 
Laplace transform of () involve integrals along the Bromwich 


contours shown in Figure 2 for o = 0 and o = o,, respectively. 
While these integrals are formally defined for infinite limits of 


integration ( jo- +00 ) , these integrals are usually approximated by 


integrals with limits ( jo- +o,,) as noted in Figure 2. 


Fourier jo Laplace 
Transform Ny, kK Transform 


JO, = Z70)f e-fpoeseeeeste eee 


Complex frequency plane 
Ss=oO+jo 


jp = -2 rf ge ference F 


Figure 2. Integration paths along the Bromwich contours for 
the inverse Fourier and Laplace transforms in the complex s- 
plane. 


The Laplace transform integrals (8) and (9) can be evaluated by 
direct integration, or as previously noted, by using the FLT which 
employs the FFT algorithm for performing the necessary 
integrations. Reference [9] provides a good review of various 
numerical quadrature methods, as well as an introduction to the 
PLT approach. 


It should be noted in passing that the FLT algorithm for 
discretely sampled data can be also described by the Chirp z- 
Transform, as developed by Rabiner in [10]. In this approach, the 
finite integration path along the Bromwich contour shown in Figure 
2 is mapped onto a contour in the complex z-plane. This technique 
is beyond the scope of this monograph, and the interested reader is 
referred to Rabiner’s work, or to refs. [11 & 12]. 


2.2. Selection of the Damping Constant 6, 


Not specified in the above discussion is what should be the 
value of the damping constant o,. At a minimum, it should be 
sufficient to ensure that (1) is satisfied for the given transient 
waveform being transformed. This implies that each of the “rogue” 
waveforms of Figure 1 could have a different value of 6,. Typically 
the value for the damping constant 6, is obtained by first examining 
the early time behavior of the computed inverse Laplace transform 
of a system response for proper causality. If a null amplitude up to 
the turn-on time of the response is not noted, then the damping 
constant should be increased. However, if 6, is increased too much, 
there will be exponentially growing noise in the late-time portion of 
the waveform, and this indicates that the value of 6, is too large. 
Thus, there is a trade-off in the accuracy of the early- and late-time 
portions of the inverse transform. 


Two suggested values for 6, based on the temporal extent of 
the waveform time window have been made in [6], and these 
provide starting points for selecting this parameter. These empirical 


estimates are o, = In (N7)/ Ta, and 6, =47/T). -, where T,,,.. is 
the maximum extent of the transient response time window, and N 


= 2”, which is the number of waveform sample points. 


In our investigations, however, we have found that these 


values of 6, provide significant errors at the end of our transient 
responses. An alternative value that we have found more useful is 


o, =K/T,_., where the value of « ranges from 3 to 7, depending 


on the nature of transient waveforms. 
2.3. Determination of a System Response 


While (8) and (9) provide the transform relationships between 
a transient function and its spectrum, it remains to define how the 
waveform or spectrum would be modified as it passes through a 
linear system transfer function [13]. Figure 3 shows a block diagram 
of such a system which has a complex frequency domain transfer 
function H(s) that relates the output spectrum G(s) to the input 
spectrum F(s) as 


G(s) = H(s)F(s). (10) 
In the time domain, this relationship is equivalent to the 


convolution function between the transient excitation f(/ and the 
impulse response of the system /(A) as 


g(t) =| hr) f(t—)dr 


= [ron —t)dt (11) 
= f(N*h) 


In this expression, the impulse response of the system /(A) is the 
inverse Laplace transform of H(s), which is given by (9). 


Figure 3. Representation of a linear system with transfer 
function H(s), excitation spectrum F(s) and _ response 
spectrum G(s). 

A unit-amplitude transfer function H(s) = 1 simply replicates 
the excitation function f(/ as the response function g(/. The next 
simplest system transfer function is one that provides a 
distortionless time delay by the amount 4 between the input and 
output waveforms, as g(/) = f(¢ — ¢). In the Fourier transform 
domain where o, = 0, this is described by the transfer function 
H(@)= exp(- jot, ).In the Laplace domain this time-shift transfer 


function is H(s)=exp (—st, ) , where s is the same complex 
frequency used in the excitation and response transforms. 

For system transfer functions that can be defined analytically in 
the Fourier domain, extending these into the complex frequency 


domain is simply done by replacing j@ by s, as noted in the time- 
shift example above. Sometimes, however, the Fourier domain 


transfer function H(@) is not represented analytically, but it consists 


of numerical values that have been obtained as a function of the 
real (measurable) frequency. This amounts to a Fourier domain 
transfer function defined along the j@ axis in the s-plane that must 
be analytically or numerically continued into the complex frequency 
plane. 


The determination of the Laplace domain transfer function 
from the Fourier domain transfer function can be done numerically 
in two different ways. The first method is to recognize that the 
Laplace spectrum given by (8) is essentially the Laplace transform 
of the product of two transient functions, (4 and exp(-sd). 
Denoting the inverse Fourier transform of /(4 as H(@), and noting 
that the corresponding Fourier spectrum of the exponential 
function is 1/(j@+0,), the convolution theorem [3] permits us to 
write the Laplace spectrum H(s) as 


H(s)=H(a)* 


jotro 


oO 


= [ (ag, or (12) 
LO @-O+6 


oO 


= | H(@-2)——aé 
es I(S)+o, 

This procedure of (12) is a convolution of two frequency 
domain spectra, and for spectral responses that have many data 
points, this calculation canl take a significant amount of time. An 
alternate and significantly faster approach for determining H(s) 
from H(@) is to use the inverse Fourier transform operator ¥' of 
(3) to determine the impulse response /(4, and then to take the 
direct Laplace transform £ of (8). This procedure is shown 
symbolically in the following sequence: 


H(s)=2£| ¥'(H(o))| (13) 


Of course, in performing this calculation, the same value of the 
damping parameter o, used for determining the Laplace transform 


of f(4) should be used in (13). 
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Chapter 3. Use of the Fast 
Laplace Transform 


In this section we provide several examples of the use of the 
PLT on waveforms and responses of interest in the EMC field. 


3.1. The FLT Applied to Simple Analytical Waveforms 


The difficult waveforms of Figure 1 provide good examples 
for the use of the FLT. In this case, we assume that the waveforms 
are passed through a system with a distortionless delay function 
H(s)=exp(-st,) with ¢, = 20 us. By applying the FLT with an 
assumed damping parameter of o, =1.524x10°to the four 
waveforms, multiplying each by H(s) and then taking the inverse 
PLT, we obtain the time-shifted waveforms shown in Figure 4. 
These waveforms are well behaved for both early and late times, 
and they are clearly a 20 ws shifted replication of the original 
waveforms. These waveforms provide absolutely no problem for 
the FLT processing. 


To contrast the FLT results with those provided by the FFT 
operating on the same waveforms in Figure 1, the FLT was re-run 
with 6, = 0, which as we have seen, provides just the FFT. The 
resulting waveforms are shown in Figure 5. All of these waveforms 
exhibit the usual fold-through effect at early time, where the non- 
zeto portions of the late-time waveforms appear in early time, when 
the responses should be identically zero. 
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Lightning oo = 1.524 10° 
channel 
current 


Lightning-produced 
E-field 


f(t) 


Sinusoidal signal 


0 20 40 60 80 100 
Time (us) 
Figure 4. Transient responses of the waveforms of Figure 1 


that have been obtained from a 20 Us time-shift operator in the 


Laplace spectral domain (o, = 1.524 x 10°), using the FLT 
procedure. 


Lightning channel 


Pel current ao= 0 


Lightning-produced 


f(t) 


Sinusoidal signal 


0 20 40 60 80 100 
Time (ps) 


Figure 5. Transient responses of the waveforms of Figure 1 
that have been obtained from a 20 Us time-shift operator in the 
Fourier spectral domain (6, = 0), using the FTT. 


It is also interesting to examine the behavior of the time- 
shifted waveforms using the FLT when a larger damping constant is 
used. Figure 6 shows the waveforms obtained for a damping 
constant of 6, = 4.998 x 10°, which is significantly larger than the 
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suggested values of Section 2.2. We note that the early time 
behaviors of the shifted waveforms are fine, but at late times there 
is noise in the response that grows exponentially. This is due to the 
fact that such a large exponential attenuation of the original time 
function f(4 in the Laplace transform of (5) has reduced the late 
time portion of the waveform to essentially the numerical noise 
level in the computation. Upon taking the inverse transform and 
multiplying by the growing exponential term of (7), the noise is 
amplified and it washes out the desired waveform. 


400 


Sp = 4.998 x 10° 


300 


Time (us) 
Figure 6. Time shifted responses of the waveforms of Figure 


1 for a large damping constant (o, = 4.998 x 10°) using the 
FLT. 


This divergent behavior of the inverse Laplace transform has 
been examined in [14], where it is stated that the inversion of the 
Laplace transform is a good example of an exponentially ill-posed 
problem. In [14] the authors also use the FFT to compute the direct 
and inverse Laplace transforms of transient functions, and they 
discusses the effects that noise and filtering have on the inverse 
transform. They seem rather pessimistic, however, about the use of 
the Laplace transform, stating that “Our results give cogent reasons 
for the general sense of dread most mathematicians feel about 
inverting the Laplace transform.” 


Notwithstanding the cautionary remarks in [14], we have noted 
that the FLT algorithm and its inverse work very well for practical 
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waveforms encountered in the EMC area. Several more examples 
will be provided in this section. 


3.2. Application of the FLT to EMC Problems 


The previous section has illustrated the use of the FLT for 
transient waveforms modified by a very simple system transfer 
function — a time shift operator. In this section we examine the 
effects of passing the signal through more complicated transfer 
functions that arise from EM field coupling to a highly resonant 
transmission line and EM field interaction with a buried facility that 
is characterized by a measured continuous wave (CW) transfer 
function. Finally we conclude with a circuit example that illustrates 
the use of the FLT on a problem with specified initial conditions. 


3.2.1. A Highly Resonant Transmission Line 


The use of FFT methods for determining the transient 
responses for the load voltages and currents on the two-wite 
transmission line of Figure 7 has been described in detail in [15]. 
The basis of this analysis is the Baum-Liu-Tesche (BLT) equation, 
which was developed in the time-harmonic (frequency) domain. By 
substituting s for the variable j@ in the BLT equation, we can 
express the load voltages 1’, and l”, at each end of a lossless 
transmission line in matrix form as 


inc 


d (1 = eri (l-cos(y ie ) 


Me 7 i 4 Pp, 0 oe ele I 5 
J : sL/c inc : 
V,(s) 0 1+ p, | |e —P> = Ed eitle (1 _ @l(l-cos(y We ) 


2 

(14) 

This equation is equivalent to (7.40) in [15], with the 

simplification that the propagation constant of the transmission line 

is equal to the free-space propagation constant. In this equation, /, 

and p, are the voltage reflection coefficients at each end of the line, 

L is the line length, d is the wire separation and E” is the incident 

E-field strength. As shown in Figure 7 this incident field arrives at 

an angle y relative to the line orientation. A similar BLT equation 

can be developed for the load currents, but the details are not 
presented here. 
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Ci Radius a 


Figure 7. A two-wire transmission line illuminated by an 
incident plane wave EM field. 


As an example of the use of the FLT for this problem, we 
consider a transmission line having the following dimensions: L = 1 
m, d = 1 cm and conductor radius a = 1 mm. The line separation 
and conductor radius define the characteristic impedance of the line 
as Z,=120 In(d/a) = 276.3 Q. This impedance is used with the 


load impedances to determine the voltage reflection coefficients, as 
p=(Z-Z)/(Z+Z). 

The excitation of this transmission line is provided by the 100 
V/m exponential waveform shown in Figure 8. This waveform has 
a 1 us decay time and a waveform start time of 0.1 Us. For this 
example, the angle of incidence of the excitation field is y= 45°. 


In this example, we calculate the voltage at load #2 and the 
current at load #1 for a highly resonant line. This is difficult to do 
using conventional FFT methods. To obtain such a resonant line, 
we will use a short circuit for the left load (Z, = 0) and an open 
circuit for the right load (Z, = 0). 
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Figure 8. The incident E-field excitation waveform for the 
transmission line. 


Reference [15] discusses how transient responses of highly 
resonant structures like this example can be solved in the spectral 
domain by expanding the inverse matrix in (14) in a suitable manner 
and then computing a modified spectrum using only a few terms of 
the expansion. The FFT of this modified spectrum provides a 
waveform that is correct for several oscillations, but then goes to 
zero for later times. This property of the waveform permits the use 
of the FFT in the analysis. 


If this spectral modification procedure is not used and the 
voltage or current spectra are calculated along the j@ axis of Figure 
2, there are periodic singularities that occur in the spectra, due to 
the resonances on the transmission line. These are shown in the 
spectral plot of Figure 9. If the FFT is applied to these spectra, 
incorrect waveforms like that shown in Figure 10 result. Note that 
the amplitude of this waveform is unreasonably large and the 
waveform is non-causal, as there is a response prior to the turn-on 
time of the excitation at 0.1 Ls. 
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Figure 9. Plots of the voltage and current transfer function 
magnitudes | V,/E*| and | I,/E”*| along the j@ axis (0, = 0). 
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Figure 10. The transient response for V, as computed from 
the voltage transfer function of Figure 9 using the FFT. 
(Incorrect) 


The waveform is physically incorrect for two reasons. The first 


is because we have neglected possible loss in the transmission line, 
the line resonances are infinite. In reality, loss will provide a finite Q 
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to these singularities, and the response waveform in Figure 10 will 
eventually die out. The second reason is due to the fact that we are 
sampling the FFT spectrum only at discrete points. If insufficient 
points are taken, the spectral peaks are not well-sampled and the 
computed transient response is erroneous. 


Thus, to calculate this response accurately using an FFT, loss 
must be added to the problem, and more points added to the 
waveform and its spectral response. Changing the Z, load to 1 Q 
will provide a suitable loss to the structure, and increasing the time 
window of the exponential excitation waveform to 6 us will 
increase it sufficiently so that the response of the line has died out. 
Then, increasing the number of points in the waveform from 4086 
to 32,768 will provide for better sampling of the resonances. 


Figure 11 shows the resulting early-time response for the 
voltage 1’, with these modifications to the problem parameters, for 
the FFT processing. Aside from the very slight decay in the 
oscillations on the line due to the added loss, the response is a 
reasonable approximation for this highly resonant line. 


V2 (Volts) 0 


0 0.05 0.1 0.15 0.2 
Time (us) 


Figure 11. Plot of the transient voltage response V> using the 


FFT by adding a small amount of loss, extending the 
waveform window and increasing the number of points in the 
waveform and spectrum. 
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Alternatively, the FLT can be applied to the computed Laplace 
spectrum for a suitable value of o, = 3.612 x 10° (and with the 
smaller sampling density of only 4096 points in the transient 
waveform). Figure 12 plots the Laplace spectral magnitudes for the 
voltage and current transfer functions along the Bromwich contour 


6, +j@. Note that in Figure 12 the abscissa variable is the frequency 
f= ©/2n, which is derived from the imaginary part of the complex 
frequency s. Figure 13 presents the computed transient voltage VV, 
for both early and late times. 

The short circuit current I; has been calculated in a similar 


manner using the BLT current equation, and its transient response 
using the FLT is shown in Figure 14. 
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Figure 12. Plots of the voltage and current transfer function 
magnitudes |V,/E| and |I,/E”| along the Bromwich 
contour 6, +j@. 
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Figure 13. Plots of the early time (left) and late time (right) 
waveforms for the voltage Vz on the lossless transmission line 


using the FLT. 
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Figure 14. Plots of the early time current waveform I, using 
the FLT. 


It should be pointed out that this FLT analysis approach also 
can be applied to transmission lines having loss, with no undue 
difficulties. This is unlike the more conventional analysis of Chang 
and Kang [16]. 
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3.2.2. Use of the FLT with a Measured CW Transfer 
Function 


Measured CW (time harmonic) transfer functions are often 
used to characterize a linear time invariant system. Such transfer 
functions have both a magnitude and phase, and if these functions 
are known over a sufficiently large frequency range, the transient 
response of the system to an arbitrary transient excitation can be 
computed. 


One example of such a transfer function has been described in 
[17], where radiated field measurements were made on an 
underground bunker. The facility was illuminated by a wide-band 
antenna, and measurements of the transfer function between a 
component of the internal magnetic field and the excitation voltage 
of the antenna were made. A network analyzer was used to measure 
both the magnitude and phase of the transfer function. 


Figure 15 shows the magnitude of a typical transfer function 
measured from 1 MHz to about 1 GHz. The corresponding phase 
for this transfer function was also measured and is used in the 
calculations, but it is not displayed. As noted earlier, this transfer 
function is measured along the j@ axis of Figure 2 and for use with 
the FLT. It must be analytically continued onto the Bromwich 
contour 6, +j@. This will be done using the procedure shown in 


(13). 


Figure 16 shows the computed impulse response /(4, which is 
computed by an inverse FFT applied to the spectrum of the 
measured transfer function of Figure 15. For accurate processing, 
such transients should always be checked for causality and for 
proper behavior at late times. This particular transient is well 
behaved — indicating that the measured CW data are sufficiently 
accurate. 


For the second step in evaluating (13), the FLT is applied to 
A with an assumed damping constant 6, = 6.62 x 10° This 
provides the value of the transfer function along the Bromwich 


contour 6, +j@. Figure 17 provides an overlay plot of the measured 
and analytically extended FLT spectra. 
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Figure 15. A typical measured transfer function magnitude for 
the internal magnetic field in a buried facility [17]. 
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Figure 16. Inverse FFT of the measured complex CW transfer 
function H(q@) of Figure 15, yielding the system impulse 
response h(t). 
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Figure 17. Overlay plot of the measured transfer function 
magnitude with the analytically continued transfer function on 
the Bromwich contour 6, +ja. 


To illustrate the FLT processing of this transfer function, we 
assume that the antenna is excited by a unit-amplitude time-shifted 
sinusoidal waveform 


f\=e™ sinQQz f,(t-t,)@(t-1,) volts 


with f = 100 MHz, a = 1 x 10” (s’) and # = 0.05 ps. Figure 18 
plots this waveform and its computed FFT and FLT spectra. 


To determine the internal transient response for the magnetic 
field for this particular antenna excitation, we could convolve the 
waveforms in Figure 16 and Figure 18, but this would take a 
significant amount of computation time. Alternatively, we can 
multiply either the FLT or the FFT spectra of the excitation by the 
appropriate transfer functions of Figure 17 and take the inverse 
PLT or FFT to obtain the transient responses. 
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Figure 18. The assumed transient voltage source of the 
radiating antenna (top) and the corresponding FLT and FFT 
spectral magnitudes (bottom). (0, = 6.62x10° for the FLT 
spectrum). 


Figure 19 presents the transient response resulting from the 
FLT processing using 6, = 6.62 X 10°. The late time response does 
not show the exponential growth of the noise, so the damping 
parameter is appropriate. There is, however, some very slight non- 
causal behavior in the early time response, which is attributed to the 
fact that the measured spectrum had to be interpolated to 
correspond to the frequencies of the FLT spectrum of the 
excitation. This interpolation of the measured data introduces small 
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errors in the measurements that give rise to this anomalous, but 
small effect. 


In examining the turn-on time of the waveform in Figure 19, 
we note that the antenna excitation waveform turns on at 0.05 Ls, 
and the response in the facility has an additional time delay of about 
0.07 ws. This provides a total time delay of 0.12 us, which is 
indicated in the figure. 


440-7 
5x10 ® 
g(t) 0 
~sao~® 
—1x1077 
0.05 0.1 0.15 0.2 
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4x1077 


5x10 


Late Time (us) 


Figure 19. Early- and late-time calculated transient responses 
of the internal magnetic field for the antenna excitation of 
Figure 18 using the FLT processing. 
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Figure 20 shows the resulting response waveform for the FFT 
processing of the spectra. It is clear that the response waveform is 
not causal. Furthermore the time from the transition of the very 
eatly non-causal sinusoidal waveform to the apparent system 
response waveform is not correct value of 0.12 us. This FFT 
calculated response is clearly incorrect, and it shows the difficulties 
encountered in trying to use an FFT process in this case. 
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Figure 20. Early- and late-time calculated transient responses 
of the internal magnetic field for the antenna excitation of 
Figure 18 using the FFT. 
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3.2.3. Use of the FLT on Problems with Initial Conditions 


A very useful attribute of the Laplace transform is its ability to 
include initial conditions in the analysis. This benefit is also found if 
the transforms are performed numerically using the FLT. To 
illustrate this, we consider a Marx pulser model with a peaking 
citcuit, as shown in Figure 21. The use of this circuit for producing 
fast-rising pulses has been discussed by Taylor and Giri [18] and we 
will illustrate the use of the FLT for conducting an analysis of this 
generator. 


The circuit consists of a Marx pulser, having circuit elements 
R, L and C, and a peaker capacitance Cy as a load. The secondary 
peaking circuit contains a load inductance L, (which is usually a 
parasitic element) and a load impedance Z;. Initially, the 
capacitance C is charged to a voltage (0+) = — V’, with the other 


capacitor and inductors having zero initial voltage and current. At ¢ 
= ¢,, the Marx switch is closed and the loop cutrent 7;(4 begins to 
1 P 1 8 


flow in the Marx section’. The current 7/(4 charges the peaker 
1 g Pp 


capacitance and when the charge is sufficient, the peaker switch is 
closed at time ¢ = ¢. This results in a very fast rising current in the 


peaker circuit and a transient excitation of Z,. 


For this example we assume the following parameters: L = 3 


WH, L, = 1nH, C= 5 nF, C,= 1 nF and R= 0.1 QQ. 


Two different values for the peaker load Z; are considered: 1) 
a fixed load of 50 Q, and 2) a load provided by the input impedance 
of a lossless 50 © transmission line with a length of 60 meters and 
loaded by a 20 Q resistance at the end. 

The switching times for this circuit are assumed to be ¢, = 0.1 


us, and fp = 0.2 ps. These parameters are not optimized for a 


particular pulser design, but are used only for illustrative purposes. 


4 In this development, we will use lower case variables like »(/ and i(t) to denote transient 
responses, and upper case variables V(s) and I(s) to denote the Laplace transform 
responses. 
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Figure 21. A Marx generator and peaking circuit used to 
illustrate the use of initial conditions in the FLT. 


The Laplace transform analysis circuits having initial 
conditions is standard material in circuit analysis texts. In [13] for 
example, it is shown that if L (fA) = F(s) represents the Laplace 
transform operation of (8), the transform of the derivative of f{(/) is 
given by 


[4] -sF0)- fO,), (15) 


where f(0;) is the initial value of the transient response as ¢ + 0, as 
approached from ¢> 0. 

Applying (15) to an inductor and capacitor, which are 
described by the relationships v(t) =L di/dt and i(t)=C dv/dt, 


respectively, the l’-I relationships for these elements in the Laplace 
domain are 


V(s)=sLI(s)—Li(O,) (for the inductor) and (16a) 
I(s)=sCV(s)—Cv(0,) (for the capacitor). (16b) 


After the Marx switch is closed, and prior to the peaker switch 
being closed at ¢ = ft, Kirchhoff’s voltage law (KVL) may be 


applied to the current loop in the Marx section to determine the 
current I7. As this requires the sum of all of the voltages across the 


citcuit elements in the loop be zero, it is convenient to invert (16b) 
to obtain the capacitor voltage in terms of the current, which is 
I(s) vO 
V(s)= OY a 2 (for the capacitor). (17) 


Thus, the KVL for the Marx section becomes 
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0 
ea | reo) : (sLI,(s) - Li,(0,)) RI(s) t its) t onl J =0 
sC s sc, s 

(18) 
Since the initial current 7;(0,) and peaker capacitance voltage 


V+) ate both zero, and since the initial Marx capacitor voltage 


is v(0+) =— V,, the loop current in (18) can be expressed as 
Ke: 
1(s)= : | (19) 
SL ee + es +R 
s\C C, 


Equation (19) represents spectrum for the waveform that starts 
at ¢ = 0. The current spectrum can be modified by the time shift 
Laplace operator exp(—st,) to account for the Marx switching time 


at f= ry. 


The usual way of taking the inverse transform of (19) to obtain 
the transient current in the Marx circuit is to expand the equation 
into a partial fraction expansion and then take an analytical inverse 
transform of each term. Alternatively, we can obtain a numerical 
solution by using the FLT on the Laplace spectrum of the current. 


Without the peaker switching, the transient current 7;(4 will 
oscillate for a long time until it is eventually damped out by the 
resistance R. When the peaker switch is activated, however, a 
second current loop is introduced into the circuit and this must be 
included in the circuit analysis. This is done by developing two 
coupled mesh equations for the spectral currents I7(s) and p(s). 


These equations involve the current through L and the 
voltages across C and C, at time fy, as initial conditions. 
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The two mesh equations can be expressed in matrix form as 


lj 1 1 1 : 1 
aofbed or ~3C,, = Li()——(Velt) + Vep(t)) 
er V-, (t,) 
ase sL, +——+Z, : ane 
sC sC Ss 


(20) 
where the right hand vector is the forcing function, and the 2x2 
matrix to the left is the system impedance matrix. This equation can 
be inverted either numerically or analytically to permit the 
evaluation of the currents I; and Ip. 


The resulting transient currents flowing in the pulser circuit 
were determined by first defining a maximum time window desired 
for the responses, which was chosen to be 20 Us. With a choice of 
2'3 = 8192 sample points, a sampling interval of 2.44 ns is obtained. 
The value for the damping parameter of the Laplace transform in 
(5) was chosen as 6, =K/T_,,, with « = 4. This gives o, = 0.2 x 
109, 

The transient response for 7,(4) was computed from the inverse 


PLT of (19). Figure 22 shows this current starting at 0.1 us and 
continuing as a slowly damped sine. The inverse FLT of I,(s) from 
(20) was computed and shifted to begin at 0.2 Us, just after the 
peaking circuit switches. This response is the dotted line in Figure 
22, and it shows how the Marx current changes when the switching 
occurs. The other solid curve in this figure represents the current 
flowing in the peaker circuit, 7(4. This response is seen to have a 
very fast rise time, which is controlled by the load resistance and the 
peaker inductance L,. For this calculation, the peaker circuit load 


was taken to be the 50 © resistive element. 
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Figure 22. Transient response of the Marx generator current 
i,t) assuming no peaker switching, and the currents i,(t) and 


i,t) for peaker switching at 0.2 ps. (Peaker load is 50 ©.) 


Once the current spectral responses are determined, the 
voltages across the capacitors can be determined using (17). Note 
that this requires including the initial conditions of the capacitor 
voltages in the inverse FLT. Figure 23 illustrates the transient 
voltages Vc and Vop as computed using the FLT. The solid curves 
represent the responses if no switching of the peaker circuit occurs, 
while the dotted curves show the behavior of the voltages after the 
switching at ¢ = tp. 


While the previous results could have been obtained 
analytically using the conventional partial fraction expansion 
method, it is more difficult to do this if the load impedance Z7 is a 
measured function, as in the previous example, or if it is a non- 
rational function of frequency. In these cases, the use of the FLT 
makes the analysis significantly easier. 
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Figure 23. Plots of the transient voltages across C and Cp» with 
and without the switching of the peaker circuit at t = ty 


(Peaker load is 50 Q.) 


As an example, we consider the case of a section of a 60 meter, 
50 © transmission line connected to the pulser instead of the 50 Q 
resistance. The input impedance of this line, when terminated in the 
20 Q load, can be calculated from the transmission line models 
given in [15]. This input impedance is shown by the solid curve in 
Figure 24. The analytical continuation of this impedance function 
onto the 6, = 0.2 x 10° contour in the complex s-plane can be 
done as in the previous example through the use of (13), and this is 
also shown in the figure. 


Figure 25 presents the impulse response of the input 
impedance of the transmission line. This clearly exhibits the 
multiple reflections occulting on the transmission line. 
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Figure 24. Plot of the input impedance magnitude |Z;,(j@) | 
and the Laplace spectral magnitude | Z;,,(s)| for the 60 meter 


transmission line load. 
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Figure 25. The impulse response of the impedance function, 
as determined by an inverse FFT of the computed function 


Z in(@)- 
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Figure 26 plots the Laplace spectra for the voltage across the 
load element Z; as a function of frequency for both the 


transmission line load and the constant 50 Q load. The presence of 
the oscillations on the transmission line is evident. 
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Figure 26. Plot of the load voltage Laplace spectra for the 
transmission line load and the 50 Q load. 


Figure 27 shows the corresponding transient responses for the 
load currents, which have been time shifted to start at time 7). It is 
seen that for a time ¢ = 2 X 60/c = 0.4 Us after the switching at ¢) 
(or at 0.6 ms absolute time) the waveform appears like that across a 
constant 50 © load. However, after the end reflection arrives back 


to the pulser circuit, the response deviates from the constant load 
value. 
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Figure 27. Plot of the transient load voltage across Z7, for the 


transmission line load and the 50 Q load. 
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Chapter 4. Summary 


This monograph has reviewed the use of the fast Laplace 
transform, as implemented using the fast Fourier transform. This 
method has been described earlier by investigators in the electrical 
power community, but it does not seem to be widely used in the 
EMC community. 


We have illustrated the use of the FLT for time shifting several 
waveforms that pose difficulties when used with the FFT. This time 
shift procedure is performed in the spectral domain by multiplying 
the Laplace spectrum by the time shift spectral operator, and then 
inverting back into the time domain with the inverse fast Laplace 
transform. The use of this method on a highly resonant 
transmission line and on measured CW transfer function data also 
has been discussed. 


Mathematicians have described the inverse Laplace transform 
as being “all conditioned” and refer to this as a “bad truth” about 
the procedure. In our work, we have observed that this ill 
conditioned nature of the Laplace inversion occurs in the form of 
exponential noise in the late time portion of the inverted transient 
response. The reason for this ill conditioning is similar to trying to 
find the limit of sin(x)/x numerically. Eventually this limit becomes 
noise dominated and the actual limit of unity is never found. 


In the inverse Laplace transform case, almost all computed 
time domain functions will become unstable if the time is allowed 
to run long enough. Hence, this method is “bad” because it does 
not replicate the actual waveform in late times. However, this 
procedure does not limit the utility of the method for early and 
intermediate times, as noted in the examples presented here. 


The key to obtaining useful Laplace transform results is 
determining a suitable value for the damping constant, o,. As noted 
here, this parameter must be positive and its “optimal” value cannot 
be defined explicitly. Several suggested values for 6, have been 
given for waveforms based on their time duration or their sampling 
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rate. But ultimately, it is up to the user of this procedure to decide 
upon a suitable value. 


Different values of 6, will provide different responses of the 
inverse transform at late times, and perhaps this ambiguity in the 
inversion process is what upsets the mathematicians. 
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absolutely integrable, 1 

angle of incidence, 15 

bilateral Laplace transform, 6 

BLT equation, 14 
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causality, 21 

Chirp z- Transform, 7 

complex frequency s, 3 

computed impulse response, 
21 

convolution, 9 

Cooley and Tukey, 1 

damping constant, 8 

Dirac delta function, 2 

discrete Fourier transform, 1 

exponential waveform, 15 

fast Fourier transform, 1 

fast Laplace transform, 4 

PLT. See fast Laplace 
transform 

FLT spectrum, 24 

Fourier transform, 5 

Gauss, 1 

highly resonant structure, 16 

ill-behaved waveforms, 3 

impedance matrix, 30 

impulse response, 9 

incident E-field, 14, 16 

initial conditions, 14, 27, 28, 
29, 31 

input impedance, 32 


inverse Fourier transform, 6 

Laplace transform pait, 6 

late time response, 24 

lightning return stroke 
current, 2 

Marx generator, 28, 31 

Marx pulser, 27 

measured CW transfer 
function, 21 

mesh equations, 29, 30 

modified Fourier transform, 4 

non-causal, 16, 24, 26 

numerical Laplace transform, 
3 

peaker switch, 27, 28, 29 

plane wave, 15 

reflection coefficients, 14, 15 

sinusoidal signal, 2 

sinusoidal waveform, 1, 23, 
26 

spectral modification 
procedure, 16 

system response, 8 

time shift, 9 

transfer function, 4, 8, 9, 10, 
14, 17, 19, 21, 22, 23 

transmission line, 14, 40 

transmission line oscillations, 
34 

turn-on time, 8, 16, 25 

underground bunker, 21 

wide-band antenna, 21 
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This monograph reviews the use of the Laplace 
transform as implemented numerically using the 
fast Fourier transform. This method has been 
described earlier by investigators in the electrical 
power community, but it does not seem to be 
widely used in the electromagnetic compatibility 
(EMC) area. The goal in developing this 
monograph is to bring this computational method 
to the attention of the workers in this community 
by providing several examples and comments on 
its use for practical problems. 
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